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Abstract. In this article, wc continue our mathematical study of organic solar cells 
(OSCs) and propose a two-scale (micro- and macro-scale) model of hctcrojunction 
OSCs with interface geometries characterized by an arbitrarily complex morphology. 
The microscale model consists of a system of partial and ordinary differential equa- 
tions in an heterogeneous domain, that provides a full description of excitation/- 
transport phenomena occurring in the bulk regions and dissociation/recombination 
processes occurring in a thin material slab across the interface. The macroscalc 
model is obtained by a micro-to-macro scale transition that consists of averaging 
the mass balance equations in the normal direction across the interface thickness, 
giving rise to nonlinear transmission conditions that arc parametrized by the in- 
tcrfacial width. These conditions account in a lumped manner for the volumetric 
dissociation/recombination phenomena occurring in the thin slab and depend lo- 
cally on the electric field magnitude and orientation. Using the macroscale model 
in two spatial dimensions, device structures with complex interface morphologies, 
for which existing data arc available, arc numerically investigated showing that, if 
the electric field orientation relative to the interface is taken into due account, the 
device performance is determined not only by the total interface length but also by 
its shape. 

Keywords: Organic solar cell; nonlinear reaction-diffusion system with electrostatic 
convection; scale transition; multiscale analysis; numerical simulation; finite element 
method. 



1. Introduction and Motivation 

Research on photovoltaic energy conversion has recently received great impulse due to 
the growing demand for low carbon dioxide emission energy sources. In particular, the 
high manufacturing cost of crystalline silicon and the latest advancements on semicon- 
ducting polymer design and synthesis in recent years have directed the attention of the 
scientific community towards Organic Solar Cells (OSCs), i.e. solar cells based on organic 
materials 12, 20 , 21, 32, 33, 38, 42 , especially because of the very limited thermal budget 
required for the production of such materials and of their amenability to be deposited 
on large areas, which is fundamental in light harvesting applications. One of the main 
peculiarities of OSCs is that most physical phenomena that are critical for charge pho- 
togeneration occur at the interface between the two materials that constitute the active 
layer of such devices. In order to increase cell efficiencies, currently of the order of about 
10% [22] , the optimization of the morphology of such interface is considered by device de- 
signers to be an issue at least as important as the optimization of the donor and acceptor 
optoelectronic characteristics [40] . 

For this reason, in this article we continue our mathematical study of organic photo- 
voltaic device models started in [IT], and we focus on the accurate and computationally 
efficient modeling of the main dissociation/recombination processes occurring in a thin 
material slab across the material interface and evaluating their impact on device pho- 
toconversion performance. With this aim, we consider a two-scale approach to OSC 
simulation that is intermediate between a continuum model [6] and a full microscopic 
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model [301 ITT] , and represents an extension to the case of arbitrary interface geometries 
of the one-dimensional model for bilayer OSC devices proposed in [?]. 

The approach is based on the introduction of two distinct levels of description of the 
physical system at hand, a micro and a macro scale, and of two corresponding mathemati- 
cal models based on classical mass balance conservation laws. At the microscale, a system 
of PDEs in a heterogeneous domain provides a full description of the excitation/transport 
phenomena occurring in the bulk regions and of the dissociation/recombination processes 
occurring in a thin material slab across the interface. The numerical treatment of the 
microscale model presents several difficulties related to the wide difference in size between 
the bulk regions and the interfacial width H . As a matter of fact, as polaron dissociation 
is assumed to occur in the first layer of polymer chains on either side of the interface, the 
length scale H can be taken to be that of the average separation between polymer chains 
which is typically more than two orders of magnitude smaller than the size of the bulk 
regions [H 151 |4"I] . 

Therefore, to obtain a computationally efficient model, we carry out a micro-to-macro 
scale transition that somewhat resembles model-reduction techniques used for porous me- 
dia with thin fractures [3T], for reaction problems with moving reaction fronts [25] and 
for electrochemical transport across biological membranes [35], and relies on a system- 
atic averaging of the mass balance equations in the normal direction across the interface 
thickness. The resulting macroscale model is a system of incompletely parabolic PDEs to 
describe mass transport in the materials, nonlinearly coupled with ODEs and transmission 
conditions localized at the hetero junction parametrized by the interfacial width H. These 
conditions account in a lumped manner for the volumetric dissociation/recombination 
phenomena occurring in the interfacial thin slab. The fact that in the macroscale model 
the interface is reduced to a zero width surface is further exploited to account for the 
local dependence of the polaron dissociation rate on the electric field orientation, which 
is the main advantage -together with the computational cost reduction- of our approach, 
as compared to previous multi-dimensional models [51 1441 [55] , 

An outline of the article is as follows. In Sect. [2] we illustrate the sequence of physical 
phenomena that lead from photon absorbtion to current harvesting in an OSC. Sect. [3] 
is devoted to characterizing the mathematical model of an OSC. In Sect. |3.l] we describe 
the geometrical heterogeneous structure of the device, while in Sect. |3. 2] we introduce the 
basic modeling assumptions on the dependent variables of the problem. Then, in Sect. |3.3] 
we present the microscale PDE/ODE model system with the initial and boundary condi- 
tions, while in Sect. |3. 4] we describe in detail the scale transition procedure that leads from 
the microscale model to the macroscale equation system. We complete the mathematical 
picture of a bilayer OSC by illustrating in Sect. |3.5| a novel model that we have devised 
for including the dependence of the polaron dissociation rate constant on the local electric 
field and on the morphology of the material interface. In Sect. [4] we briefly comment 
about the numerical methods used for discretizing the macroscale model, while Sect. [5] is 
devoted to presenting and discussing numerical results. In particular, in Sect. |5.l] we suc- 
cessfully perform the validation of the accuracy of the macroscale model with respect to 
the microscale model through the numerical simulation of a one-dimensional OSC under 
different working conditions. Extensive simulations of two-dimensional OSC structures 
are instead reported in Sects. [5T2] to |5.4| in order to both validate the proposed macroscale 
model with respect to previously available numerical results and to analyze its effectiveness 
in the study of complex interface morphologies. Finally, in Sect. [6] we draw some conclu- 
sions and sketch possible directions for further research in the area of OSC modeling and 
simulation. 
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2. Basic Principles of Photocurrent Generation in OSCs 

In this section, we describe the basic principles of photocurrent generation in OSCs 
only to the extent strictly needed for understanding the naming conventions adopted in 
the following sections. For a more thorough introduction to the subject, we refer the 
interested reader to 20 38, 32]. The typical structure of an OSC is constituted by a thin 
film a cross-section of which is schematically represented in Fig. [I] 




Figure 1. Structure of an organic solar cell. 



The photoactive layer of the device consists of two materials, one with higher electron 
affinity (the "acceptor", for example F8BT, P3HT) and one with lower electron affinity 
(the "donor" , for example PFB, PCBM), sandwiched between two electrodes, one of which 
is transparent to allow light to enter the photoactive layer while the other is reflecting in 
order to increase the light path through the device. 

The sequence of physical phenomena that leads from photon absorption to current 
harvesting at the device contacts is represented in Fig. [2] 
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Figure 2. Flow-chart of the photoconversion mechanisms in an OSC. 



Absorption of a photon in either material produces an electron-hole pair, usually re- 
ferred to as an exciton whose binding energy is of the order of about 0.5 -j- leV. Excitons 
may diffuse through the device until they either recombine or reach the interface between 
the donor and acceptor phases. If this latter event occurs, the exciton may get trapped at 
the interface in such a way that its electron component lays in the high electron affinity 
region while the hole component lays in the low electron affinity region. Such a trapped 
excited state is referred to as a polaron pair or geminate pair 4, 44 36, 34 and has a 
lower binding energy compared to that of the exciton state, as the Coulomb attraction 
between the electron and hole is reduced by the chemical potential drop between the two 
materials. The polaron binding force may be overcome by the electric field induced by the 
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small built-in voltage between the metallic contacts thus leading to the formation of two 
independent charged particles (one electron and one hole), otherwise the polaron pair may 
return to the untrapped exciton state or recombine. Free charge carriers move by drift 
and diffusion mechanisms and, unless they are captured along their path by the coulombic 
attraction of an oppositely charged particle and recombine at the interface to form a new 
polaron pair, they eventually reach the contacts thus producing a measurable external 
current. 



3. Mathematical Model 

In this section, we propose a PDE/ODE model of photoconversion and charge trans- 
port mechanisms in an OSC. The model relies on a two-scale approach that is based on the 
introduction of two distinct levels of description of the physical system, namely, a micro 
and a macro scale, and of two corresponding mathematical equation systems based on 
classical mass balance conservation laws. The construction of the model proceeds through 
four steps. In Sect. |3.1| we describe the geometrical and heterogeneous structure of the 
device which consists of two bulk regions (the acceptor and donor phases) separated by an 



interface region of (finite) thickness 2H, while in Sect. 3.2 we introduce the basic mod- 
eling assumptions on the dependent variables of the problem. In Sect. |3.3| we introduce 
the microscale PDE/ODE model system of conservation laws that governs transport of 
the various species throughout the device, together with its initial and boundary condi- 
tions and the generation/recombination mechanisms that occur in each subdomain of the 
heterogeneous device. In Sect. |3.4[ we describe in detail the scale transition procedure 
that is applied to the microscale model in order to obtain the macroscale equation system. 
This latter system basically consists of the same equations as in the microscale model, but 
satisfied in the separate acceptor and donor phases, coupled through a set of flux trans- 
mission conditions across the material interface that synthetize in a "lumped" manner 
the dissociation and recombination mechanisms that actually occur in the thin volumet- 
ric slab of width 2H surrounding the interface itself. The resulting macroscale model 
system is a compromise between a continuum model and a full microscopic model, and 
represents a consistent mathematical rationale and generalization of the various models 
proposed in [4l 1431 144] . We conclude our mathematical picture of the OSC by illustrating 
in Sect. |3.5| a novel model of the polaron dissociation rate properly devised for including 
the dependence on the local electric field and on the morphology of the material interface. 

3.1. Geometry of the Heterogeneous Computational Domain. A schematic 3D 



picture of the OSC is illustrated in Fig. 3(a) The device structure Q is a parallelepiped 



shaped open subset of R divided into two open disjoint subregions, Q n (acceptor) and 
Q p (donor), separated by a regular oriented surface F = 9f2 n PI dQ p [14] on which, for 
each x £ V, we can define a unit normal vector fv(x) directed from Q p into S7 n - The 
top and bottom surfaces of the structure are mathematical representations of the cell 
electrodes, cathode and anode, denoted as Fc and Fa, respectively, in such a way that 



dVt n = r c urur„ and dQ,p = T.4 U T U T p (see Fig. |3(b)[ ). We also denote by v the unit 
outward normal vector over the cell boundary d£l. 

Following [U 1431 144] . it is convenient, for modeling purposes, to associate with the 



interface F the subrcgion S1h C f2 depicted in Fig. 4(a) and defined as follows. For each 
point x € F, let t m = {x + £i>>r (x) : |£| < H} be the "thickness" associated with x. Then, 
set 

(1) 0,h = |J t x = {y e n : dist(y, V) < H} . 



The subregion fin is thus a 3D thin layer of thickness 2H surrounding F which repre- 
sents the device volumetric portion where the dissociation and recombination mechanisms 
of Sect. [2] are assumed to occur. It is worth noting that the width H is, in general, an 
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Figure 3. Geometry of the cell bulk. 
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Figure 4. Geometry of the cell bulk and interface region. 
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unknown of the physical problem. As such, it depends on x and t, it may assume different 
values in the two material phases of the photoactive layer and might in principle locally 
depend on the electric field. According to the data provided in [U 1431 144] , we assume 
for simplicity H to be a constant quantity. Based on the definition Q, we can introduce 
the two portions Q' n — Q n \ Qh and Q' p = Q, p \ Qh (see Fig. 4(b)| . Consistently, we also 
introduce the boundary portions T' n and F' p and set r± = {x ± Hvr (x) : x G T}, in such 



a way that dQ,' n = V c U F + U T' n , dn' p = Fa U V. 
T H = (r„ur p )\(r^urp). Notice that, unlike V, 
as "mathematical" interfaces. 



U r' p and 8Q.H = r + ur_UT H , where 
the surfaces T_ and T+ can be regarded 



3.2. Modeling Assumptions. Let us denote by e, P, n and p the volumetric densities 
of (singlet) excitons, polaron pairs, electrons and holes, respectively, and by J e , Jp, J„ 
and J p the associated particle fluxes. Based on the physical working principles of an 
OSC illustrated in Sect.|2j on the heterogeneous geometrical decomposition of the device 
introduced in Sect. [3~T] and the extensive numerical simulations reported in [8j, we make 
the following modeling and geometrical assumptions: 

A.l: excitons can be generated at any position in the cell, so that e = e(x,t) is a 
nonnegative function over all the cell domain f2; 

A. 2: electrons (holes) are not able to penetrate the donor (acceptor) material beyond 
the interface layer Qh, so that electrons (holes) are nonnegative functions over 
£l' n U Qh (fip U Qh), and are identically equal to zero in Q' p (Q' n ); 
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A. 3: polarons are trapped and immobile in the interface region Qh, so that P is a 
nonnegative function over Qh and identically equal to zero in fl' n U Q' p ; 

A. 4: the OSC is in the "off" state at t — CP (that is, before illumination), so that 
the initial condition for all the involved densities is e(x, 0) = P(x, 0) = n(x, 0) = 
p(x, 0) = for all x G Q; 

A. 5: the geometry of the device is an infinite periodic repetition of the computa- 
tional domain of Fig. |3j so that periodic boundary conditions are enforced for all 
variables on the lateral boundary of f2. 

3.3. Microscale Model. In this section, we illustrate the microscale model we advocate 
in this work to be a mathematical representation of the functioning of an OSC. 
We take excitons to obey: 

(2a) + V • J e = Sf + S" in ft 

where we use assumptions A.l and A. 3 to define 
(2b) Sf = Q- — in ft 




in O n U ft p 



in Q 



ii- 



and 

(2c) 

Ten. 

The superscripts B and H represent the fact that the corresponding volumetric production 
terms are defined in the bulk and interface regions, respectively. The term Q denotes the 
rate at which excitons are generated by photon absorption and is henceforth assumed to 
be a nonnegative given function of time and position while r e is the exciton lifetime in 
the bulk materials. In the interface region Qh additional dissociation and recombination 
mechanisms are taken into account and T^ iss and r]k rec represent the rate constants for the 
transition of excitons to the polaron state and that of polarons back to the exciton state, 
respectively. In particular k rec denotes the total rate of polaron recombination events and 
< i) < 1 the fraction of such events which produce a singlet exciton. As excitons have 
zero net charge, their flux is driven by diffusion forces only, i.e. the flux density may be 
expressed as 

(2d) J e = -D e Ve in Q 

D e being the exciton diffusion coefficient. At the contacts we assume perfect exciton 
quenching [37] so that 

(2e) e = onr c ur A . 

Because of assumption A. 2, the following equations hold for electrons: 



(3a) 

and for holes: 




v • j„ = s£ p in n \ n' p 

in ft' 



dp + V • J p = SZ P in ft \ ft'n 



(3b) I dt 

\ p = in n' r 



where the term S n , p is defined as 



(3c) S np 



kdissP — ■yn-P m ftff 

in Q n U Q'p. 
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Notice that S^.p is identically zero in the bulk region Q' n U Q' p as electrons and holes can 
only recombine with each other, so no recombination occurs where either of the two species 
is missing. In the interface region Qh both electrons ad holes exist so the terms in S„ iP 



take into account for polaron pair dissociation with kdias rate constant (see Sect. 3.5 for 
the model) and bimolecular recombination with rate constant 7. As electrons and holes 
each bear a non-zero net charge, their flux is driven by both diffusion and electric drift 
forces [25] . therefore: 

(3d) J n = —D n Vn — /i„nE 

and 

(3e) J p = — DpVp + p p pE 

where E is the electric field while D n , (i n and D p , p p are the diffusion coefficient and 
mobility for electrons and holes, respectively. 

Because of assumption A. 2, the following boundary conditions hold at the artificial 
interfaces separating the donor and acceptor bulk phases from the thin slab region fi#: 

(3f) v v ■ 3 n = on T_ 

and 

(3g) u T ■ J p = on T + . 

At the contacts we impose the same Robin-type boundary conditions as described in [101 
HZ]: 

(3h) — K n v ■ J n + a n n = f} n on Tc 
and 

(3i) — KpV ■ J p + a p p — P P on Ta, 

where k„, k p , a n , a p fi n , j3 p are nonnegative coefficients. 



The electric field E in (3d I and ( |3e[ ) is connected to the electric potential (p by the 
quasi-static approximation 

(4a) E = -Vip in Q, 

and satisfies the Poisson equation 

(4b) V ■ (eE) = p in Q, 

where p is the space charge density in the device. Using assumption A. 2, the piecewise 
smooth definition of p turns out to be: 

{— qn in fl' n 

q(p - n) in Q. H 
+qp in Q'p, 

q denoting the quantum of charge. The electric permittivity e is equal to e r £o, £r and 
£0 being the relative material and vacuum permittivities, respectively, with e r = e r , a in 
the acceptor phase and e r = e rjC j in the donor phase, so that e may be discontinuous 
across the interface F. Dirichlet boundary conditions for the electric potential are set at 
the contacts Ya and Yc, as follows 

(4d) (p = on T c 

and 

(4e) ip = Vappi + V bi on Y A 

where Vbi = ($a — &c)/q is the built-in voltage of the cell, 3>a and $c are the contact 
metal work functions while V app i is the externally applied voltage. 



8 



C. DE FALCO 1 ' 2 , M. PORRO 1 ' 3 , R. SACCO 1 , AND M. VERRI 1 



Because of assumption A. 3, the flux Jp is identically equal to zero in all SI and for all 
t > 0, and polarons satisfy the following ODE in the interface region 

dP e 

(5a) — = V ■ynp - (k diss + k rec ) P in Q. H 

while their density is identically zero in the bulk 

(5b) P = infi^Ufip. 



3.4. Micro-to-Macro Scale Transition. The microscale model of a bilayer OSC de- 
scribed in Sect. |3.3] can be subdivided into three distinct groups of equations: 

1) : parabolic PDEs enforcing mass conservation of excitons, electrons and holes; 

2) : an ODE describing the kinetics of photogenerated polaron pairs; 

3) : an elliptic constraint enforcing Gauss theorem in differential form to be satisfied 
at each time t > throughout the whole cell domain. 

The markedly spatially heterogeneous nature of the problem may be quite impractical 
for numerical simulation, in particular when devices with complex interface morphology 
in multiple spatial dimensions are considered. For this reason, in this section we propose 
a scale transition procedure which allows us to derive a macroscale model that is more 
amenable to numerical treatment. Other examples of multiscale mathematical approaches 
that are based on the scale separation concept and scale transition can be found in [311 

HH1 GUS- 
TO construct our multiscale model of an OSC, we abandon the perspective focused at 
the nanoscopic characteristic level adopted so far, and prefer to look at the cell from a 
"larger" distance. By doing so, necessarily, we loose control of the details (i.e, we cannot 
distinguish the region Hh from the two bulk regions fi n and fl p ), but, at the same time, 
we gain the advantage of not needing to resolve the interfacial bulk region across F. The 
resulting macroscale problem is thus posed in the partitioned domain Q \ T = Q n U fl p 
(as a matter of fact, we are still able to neatly distinguish the interface separating the 
two material phases!) without including the interfacial production terms in the mass 
balance and kinetics equations 1) and 2) introduced above. 

Of course, we cannot simply limit ourselves to neglecting these latter terms, rather, we 
do need to incorporate their effects, in the macroscale model, in an alternative way. For 
this, the simplest approach to micro-to-macro scale transition consists of replacing , at 
each point of F and for each time level, with its average a?\ across the thickness of Q,h in 
the normal direction, and then, of using o~F\ as a source term for suitable flux transmission 
conditions, to be enforced on the interface F in the case of mass balance equations. In 
the case of polaron pair equation, the averaging procedure automatically transforms the 
volumetric kinetics balance within Qh into a surface kinetics balance over F. In any case, 
the scale transition results in the introduction of suitable interfacial terms that replace 
in a "lumped" manner the volumetric dissociation/generation phenomena microscopically 
occurring in Qh- Having characterized the averaging procedure for equations 1) and 2), the 
(macroscale) differential Gauss theorem 3) remains automatically (formally) unchanged 
and is expressed in terms of the (macroscale) space charge density as in Eqns. Q. 

3.4.1. Derivation of the Macroscale Equations. The macroscale model for excitons reads: 

(6a) ~ + V • J e = Sf in fi \ Y 

with 



(6b) 



Sf = Q - — in Q \ T 
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and subject to the interface conditions: 

f [i/r-JJ=ff? our 

\ [e]=0 our 

where [/] := /„ — f p denotes for any function / : SI — > R the jump of / across the interface 
r, f n and f p being the traces on V of the restrictions of / from fi„ and Q p , respectively. 
The continuity of e at the interface is a requirement consistent with the elliptic regularity 
of both microscale and macroscale problems. 

The interfacial source term is defined as 
(6d) 

rjkrecP d£ = nk rec / Pd( / e d^ ~ nk rec P e|r- 



Tdiss J J —H 7~diss J — H 



1 1 S S 



In the above relation, ejr is the (single-valued) trace of e over F, while P is the areal 
density of the bonded pairs, defined as 

(6e) P(x,t)= [ P(x + £t x ,t)d(, VsGT 

J-H 



II 



and the midpoint quadrature rule is used for approximating the third integral in (6d I. The 
macroscale model for excitons is completed by the constitutive relation (2d I for exciton 
flux density and by the perfect exciton quenching boundary conditions \2e\. 
The macroscale model for electrons reads: 

(On 
— - + V • J„ = in £l„ 

dt 
n = in Q. p 

subject to the interface/boundary condition 

(7b) w ■ J« = crZp on r. 

The interfacial source term a^ p is defined as 
(7c) 

H r-H r-H _ 

(kdissP ~ 7 n f) d£ ~ k diss \ T / P d£- -ynpd£ ~ k diss \ T P-2H j\ r n\r p\r 

-H J-H J-H 

where definition ( |6e[ | is used in the first integral while the midpoint quadrature rule is 
again used to approximate the third integral in ( |7c[ |. The macroscale model for electrons 
is completed by the constitutive relation ( 3d I for electron flux density and by the Robin- 
type boundary condition ( |3h[ ) . 

Proceeding in a completely analogous manner as done with electrons, the macroscale 
model for holes reads: 

(8a) { |+ V ' J ^° 

I p = in Cl n 

subject to the interface/boundary condition 

(8b) vr ■ J P = —<Jn, P on T. 

The macroscale model for holes is completed by the constitutive relation ( |3e| ) for hole flux 
density and by the Robin- type boundary condition pl| . 

The conditions ( 7b I and ( 8b I assume an interesting physical meaning upon introducing 
the electron and hole current densities, defined respectively as j n := —qJ„ and j p := +qJ p , 
and the total (conduction) current density j := j n + j p . Recalling that n — (p = 0) in 
fip we have: 

j n in fi n 

j p in Sip 
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from which we get 



(8c) 



K -31=0 on r, 



that expresses the property of current conservation across the interface F. 



Integration of ( 5a I across the interface thickness yields the following macroscale model 



for the areal density of polaron pairs 



(9a) 



BP H r 

— =a P onT 



where 



(9b) 



1H ~ 
cr P = e|r + 2if 7|r n|rp|r — (kdiss\r + k rec ) P- 

Tdiss 



The macroscale model for the differential Gauss theorem is expressed by the following 
Poisson problem in heterogeneous form: 



(10a) 



V • (eE) = p in fi \ T 



with 



(10b) 



P = 



-qn 
-qp 



in Q n 
in Q„ 



and subject to the interface conditions: 



(10c) 



\v T ■ eE] = on V 

M = on r. 



Two remarks are in order with system (10 1. First, we notice that the Gauss theorem in 



differential form ( 10a I looks formally identical to the corresponding microscale formula- 



tion ( 4b I , the difference between the two methodologies being in the definition of the space 



charge density p (compare (Hq) with ( 10b l). Second, the transmission conditions ( 10c 



press the physical fact that the normal component of the electric displacement vector and 
the electric potential do not experience any discontinuity at the material interface, as is 
the case of the microscale formulation. 
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3.4.2. Summary of the Macroscale Model. For sake of convenience, we summarize below 
the macroscale model of an OSC written in primal form: 

(11a) ( 14 = °. I-T ■ D e Vej = V k rec " 
e = 
t e(as,0) = 0, 



Tdi. 



infi n UO p Efl\r 

on r, 

on r c ur A , 

Vcc G O, 



(lib) 



(11c) 



(lid) 



dP _ 2H_ 
P(as,0) = : 



i + 2Hjnp — (kdiss + k rec ) P 



— V • (D n Vn — fi„n\7ip) = 

fr • (D n \7n - ^i n n\7ip) = -k diss P + 2H-ynp 
K n v ■ (D n Vn — n n nV<p) + a n n = f3 n 
n(x,0) = 0, 



: 



^ - V • (D P V P + fi pP V<fi) 

—Vr ■ (D p Vp + fippX7tfi) = —kdissP + 2H-ynp 
k p u ■ {D v Vp + LippVtp) + ctpp = p p 
p(x,0) =0, 



on r, 
Va; G r, 

in fl n 

on r, 
on r c , 
Vx g n, 

in Jlp 

on r, 

on Ta, 
Va; G fl, 



(lie) 



—V ■ (eV(/p) = — gn 
-V ■ (eV<p) = +gp 
M = [i* ■ eV^] = 
V? = 

<P = V" app ; + Vm 



in fi n 
in Q p 

on r, 
on r c , 

on V A - 



System ( 11 1 is completed by periodic boundary conditions on r„ur p , as stated in assump- 



tion A. 5. For the physical models of the coefficients in system ( 11 \ we refer to [41 1191 123] . 
except for the description of the polaron dissociation rate constant kdiss which is addressed 
in detail in Sect. |3.5| In particular, for the carrier mobilities, we neglect the effect of ener- 
getic disorder, so that they can be assumed to depend only on the electric field magnitude, 
according to the Poole-Frenkel model. As for diffusivities, in the computations of Sect. [5] 
Einstein relations 

(12) D n = (KgT / q)n n , Dp = (K B T/q)n P 

are assumed to hold, although the proposed multiscale formulation remains unchanged if 



such an assumption is removed. In ( 12 1, Kb is Boltzmann's constant and T is the absolute 



temperature. Finally, for the bimolecular recombination rate constant 7 a Langevin-type 
relation is used [4]. 

3.5. Model for the Polaron Dissociation Rate. Numerical simulations as those re- 
ported in Sect. [5] show that the polaron dissociation rate kdiss has a significant impact on 
the cell photoconversion efficiency, for this reason we devote this entire section to modeling 
the dependence of kdiss on the electric field and on the morphology of the material inter- 
face. A commonly used polaron dissociation rate model is the Braun-Onsager model [7] 
which is derived assuming the OSC bulk to be a homogeneous medium and takes into 
account only the magnitude of the electric field. In [I] the authors propose a model for 
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fcdi SS (E) tailored for bilayer devices which is derived by performing an average over the 
admissible range of the escape angle relative to the electric field direction. In this latter 
model, the electric field is assumed to be always directed orthogonally to the interface con- 
sistently with the planar geometry of the device considered therein. The authors of [44] 
apply the dissociation rate model of H] to more complex geometries by performing an 
average along the interface of the field component normal to the contacts which amounts 
to neglecting the effect of local electric field orientation. 

To construct a novel model which also takes into account this latter effect we repeat 
the derivation of [4] with two differences. The first difference is that of removing the 
assumption that the field is normal to the interface. The second difference is that of 
considering a limited range of admissible escape directions to account for the physical fact 
that polaron pairs tend to be aligned with the gradient of the electron affinity due to the 
different materials in the two device subregions. 




Figure 5. Geometrical notation of the quantities involved in polaron 
dissociation at the material interface. 

Referring to Fig. [5] for the geometrical notation, we let 

r2TT fx/2 

(13) fc dl3S (E) = k Maa (0) / dip w(9,ip) P(E-r)d6, 

Jo Jo 

where kdiss{0) is the zero-field dissociation rate constant, r is the escape direction of the 
electron part of the polaron at the point x £ T, w is a nonnegative weight representing the 
probability distribution of admissible escape directions, and such that J Q 27r dip Jq w(0,ip) dd 
1, and P is an enhancement/suppression factor given by the Poole- Frenkel formula: 



(14) 0(z) 



e~ Az z>0 



s 2 ^^ z < 0, 

having set A — (4iTE)~ 1 q 3 (Kt, T)~ 2 . The product E • r can be expressed in terms of the 
normal component E n and the tangential component E t of the electric field as 

E • r = E n cos + E t sin 6 cos ip. 

To specify an expression for w we assume an escape direction r to be admissible only when 
the angle it forms with respect to the normal unit vector v is not too large. Indicating by 
Omax the maximum admissible value for 9 and allowing all admissible values to be equally 
likely, we obtain: 

sin 6 



w(6, tp) = < 



2vr(l 



< < 



TV 

@max < < — . 
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Two limits are of particular interest, 



0+ and Or, 



it/2. 



In the first case, Eq. ( 13 I can be checked to yield 
(15) k dlss (E) = k dlss (0) 13 (E n ). 

This corresponds to assuming that all geminate pairs are exactly aligned with the interface 
normal unit vector, thus neglecting any possible variability in their orientation due to, e.g. 
interface surface roughness and/or thermal vibrations. 
In the second case, Eq. ( 13 1 becomes 

fir/2 



(16) 



,(E) = k dl3S {0) 



d<l> 



sin f 

~2^T~ 



P (E • 



)d6>, 



which, in the special case where E± = 0, coincides with Eqs. (17)-(21) of A . Notice 
that if Et 7^ the choice 9 m ax = it/2 may overestimate the effective dissociation rate 
as it corresponds to completely neglecting the alignement of the geminate pairs with the 
electron affinity gradient. This is observed to give rise to non-physical effects as shown by 
the simulations of Sect. |5.2| Therefore, for practical purposes, the quantity 9 m ax should 
be used as a fitting parameter to be calibrated on experimental data. 

Fig. |6| s hows the dissociation rate constant (normalized to kdiss(0)} computed by 
model ( |15| l (left) and ( |16[ ) (right) for several values of the angle between E and v and 
having set T — 300 K and e r = 4. We notice that the dissociation rate computed by 
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Figure 6. Comparison between models ( 15 1 and ( 16 1 for various angles 
between E and v. 



model (16 1 has a significantly smaller range of variability than predicted by model (15 1. 
A possible explanation to this difference is related to the smoothing operated by the in- 
tegral in (16 1. The higher variability of the dissociation rate translates into an higher 



sensitivity of model ( 15 1 to the inclination of the electric field with respect to the interface 



normal as will be further discussed in the numerical results section when commenting 



Fig. 14(b) A discussion of the impact of ( 15 1 and ( 16 1 on the model predictions will be 
carried out in Sect. [5J2] 



4. Numerical Approximation 

In this section we describe the numerical techniques used to solve the mathematical 
models introduced in Sects. [3~3| and [374] The full details of the discrete system of linear 
algebraic equations resulting from problem approximation are postponed to [A] As for the 
simulation of the model in the transient regime carried out in Sect. |5.l] we have adapted 
to the case at hand the numerical method described in [T7] based on Rothe's method 
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and on the use of adaptive Backward Differentiation Formulas (BDF). In the steady-state 
simulations illustrated in Sects. |5.2| |5.3| and |5.4| all partial derivatives with respect to 



time t have been dropped out in system (111 in such a way that Eq. ( lib I reduces to an 
algebraic constraint. 

The numerical strategy adopted in the present paper is basically composed of three 
steps: 

(1) Linearization 

(2) Spatial discretization 

(3) Solution of the linear algebraic system 
Step (1) 

For model linearization, we adopt a quasi-Newton approach similar to that used in [17] , 
where, in the computation of Jacobian matrix entries, the dependence of mobilities and 
polaron pair dissociation rate on the solution is neglected. 
Step (2) 

Similarly to [17] , for the spatial discretization of the sequence of linear systems of PDEs 
stemming from Step (l)we adopt the Galerkin-Finite Element Method (G-FEM) stabilized 
by means on an Exponential Fitting technique 3 , 18 : 45 , 28 in order to deal with possibly 
dominating drift terms in the continuity equations. A peculiarity of the heterojunction 



model (111 as compared to the homogenized model of |17| is the presence of non-trivial 
interface conditions at the donor-acceptor interface, which are taken care of by means of the 
sub structuring techniques described, e.g., in 1391124) which turn out to be of straightforward 
implementation in the adopted G-FEM method. 
Step (3) 

To solve the linear algebraic systems arising from problem discretization, we employ the 
Unsymmetric Multi Frontal method implemented in the UMFPACK library ,15 as on 
current hardware architectures memory constraints are not the main limiting factor and 
the use of a direct sparse solver has the advantage of being more robust than iterative 
approaches with respect to coefficient matrix conditioning 

5. Simulation Results 

In this section we carry out an extensive computational study of the micro and macroscale 
models introduced in Sects. [3~3] and [3~4] In Sect. |5.l| one-dimensional transient simulations 
under different working conditions are carried out to verify the accuracy of the macroscale 
model with respect to the microscale system. In Sects. |5.2[ |5.3| and |5.4[ computations 
are performed in steady-state conditions, in order to validate the macroscale model by 
comparison with available results in the literature. The numerical schemes of Sect.|4]have 
been implemented in Octave using the Octave-Forge package bim [TB] for matrix assembly. 



5.1. Numerical Validation of the Accuracy of the Macroscale Model. In Sects. [374 
and |3.3[ we have illustrated two models of the operating principles of bilayer OSCs at two 
increasing levels of detail, corresponding to the macro and micro scales, respectively. The 
two modeling descriptions are expected to provide a correspondingly more refined level of 
quality in the representation of the principal physical phenomena that govern the func- 
tioning of an OSC, at the price, however, of a substantial increase in implementation com- 
plexity and computational effort, especially in the case of multi-dimensional simulations 
(mesh generation of complex interface morphologies, solution of large algebraic systems 
with possibly badly-balanced matrices). The natural question that arises at this point of 
the discussion is whether the macroscale formulation of Sect. |3. 472] is capable of returning 
an output picture of the performance of a bilayer OSC with sufficient accuracy compared 
to that of the microscale formulation of Sect. |3.3] By construction of the two models, the 
extent by which the word "accuracy" is mathematically identified cannot certainly refer 
to a pointwise comparison between micro and macroscale solutions (they will certainly 
look different!), rather, it should be concerned with average Quantities that best represent 
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Parameter Symbol Numerical value 

Acceptor relative dielectric constant e r , a 2.5 

Donor relative dielectric constant e r ,d 2.5 

Built-in voltage Vu -0.6 V 

Temperature T 298 K 

Electron mobility fi n 4 • 10~ 8 m 2 V" 1 s" 1 

Hole mobility fi p 2 • 10" 8 m 2 V" 1 s _1 

Exciton diffusion coefficient D e l-10 _7 m 2 s — 1 

Exciton lifetime r e 1 • 10 -9 s 

Exciton dissociation time t<«ss l-10~ 12 s 

Polaron pair recombination rate constant k rec 1 • 10 6 s _1 

Singlet exciton recombination fraction 77 0.25 

Polaron pair dissociation rate constant (V app i = V) kdiss 1 ■ 10 7 s _1 

Polaron pair dissociation rate constant (V app i = —Vu = +0.6 V) kdiss 2 • 10 5 

Bimolecular recombination rate constant 7 1 • 10~ 
Boundary condition parameters for electrons 



Boundary condition parameters for holes 



-1 



K n 

a n 1ms -1 
Pn 0m- 2 s 



K p 

a p lms 
f3 p 0m" 



Table 1. Model parameter values used in the simulations of Sect. |5.l| 



the overall performance of the device. In this respect, the verification test we are going to 
carry out later on, is a check of the total current per unit area jtot(t) predicted by the micro 
and macro formulations, where jtot(t) = |(i(t) + d(eE)/dt) ■ v\r cont > Pcont being either of 
the two contacts Tc or Fa- The choice of jtot for model validation is due to the fact that 
the total output current density is an easily accessible quantity in experiments, and thus 
represents the most significant parameter for assessing the photoconversion performance 
of a solar cell [32] , 

For sake of computational simplicity, we consider a biplanar OSC, so that the resulting 
spatial geometrical description can be reduced to a one-dimensional model. The total 
length L ce u of the device is equal to 100 nm, with the two regions Q n and Q p occupying 
each one half of the cell. All model coefficients are assumed to be constant quantities, and 
their values are listed in Tab.[l] We start by simulating the cell turn-on transient at short 
circuit condition (V app i = V), which corresponds to computing the device response to an 
abrupt variation at t = in the photon absorption rate Q from zero to 10 25 m- 3 s _1 . The 
simulation time interval is taken wide enough for the device to reach stationary conditions. 
In Fig. |7(a)"| the relative discrepancy between the stationary values of jtot computed with 
the two methods is reported for several values of the interface width parameter H. Results 
allow us to conclude that in the physically relevant range of variation of H, the relative 
discrepancy between the micro and macroscale models remains consistently below 10% and 
that, as expected, the predictions of the two models tend to become undistinguishable as H 
tends to zero. In Fig. 7(b) we set H = 0.25 nm and we show the time evolution of the total 
current per unit area in the two biasing conditions V app i = 0V and V app i = — Vu = 0.6V 



this latter being the "flat band" voltage. Accordingly to Fig. 7(a) in both regimes the 
curves almost coincide over the whole simulation time interval. 
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Figure 7. Comparison between microscale and macroscale models. 

5.2. Model Validation through Comparison with Existing Simulation Data. In 

this section, we aim to compare the predictions of our macroscale model to those of [431144] 
and to investigate the impact of the model for kdiss proposed in Sect. |3.5] on the simulated 
device performance. We consider the same device as in [431 1441 [13] where the acceptor and 
donor materials are F8BT and PFB, respectively. The values of the model parameters are 
listed in Tabled 

The device morphology, shown in Fig. [8] is an interpenetrating rod-shaped structure 
of donor and acceptor materials with L ce ii = 150 nm, Lr = 79 nm and 

Wr — 6.25 nm. Throughout this section, we denote by y the direction between the two 
electrodes Fc and Fa- 



w R 



Figure 8. Internal morphology with rod-shaped donor-acceptor interface. 



In [431 144) an optical model has been used to determine the exciton generation term Q. 
Here, instead, we follow a simpler approach by considering Q to be constant in the entire 
device structure and equal to the value obtained averaging the result in |43l I44| . The 
choice K n = k p = corresponds to enforcing Dirichlet boundary conditions for the carrier 
densities at the contacts and amounts to neglecting the dependence of charge injection on 
the electric field and assuming an infinite recombination rate at the contacts. 

Fig. 9(a) shows the current density-voltage characteristics in the case of an exciton 
generation rate Q — 1.53 • 10 23 m -3 s _1 . The three curves correspond to the use of three 
different expressions for the polaron dissociation rate constant kdiss, identified as follows: 
(A) the model proposed in [431 I44| with E y — \T\~ f T E y dx as the driving parameter 
for polaron pair dissociation (solid line); (B) the model (dash-dotted line); (C) the 
model |T5|) (dashed line). The result computed using model (A) is in excellent agreement 
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Figure 9. Comparison of the current- voltage characteristic lines with 
two different values for the exciton generation rate. Solid line curve 



refers to model proposed in [43], dash-dotted line refers to model (161 



while dashed line refers to model ( 15 1 



with that of Fig. 7(right) in [44] despite the above mentioned modeling differences. Model 
(A) does not account for the orientation of the electric field with respect to the donor- 
acceptor interface and is expected to overestimate dissociation in the case where E • i/p — 
0. This is confirmed by the curve for model (B). As a matter of fact, in this case all 
the dissociation directions are assumed to be equally likely and the computed output 
current density before fiat-band condition occurs (V app i < 0.6 V) is smaller than predicted 
by the solid line curve. For V apP i > 0.6 V the computed current-voltage characteristic 
exhibits a non-monotonic behavior. This latter behavior is not observed in any of the 
experimental measurements we are aware of, and is most probably to be ascribed to a too 



important contribution of the tangential component of the electric field E t that leads (16 1 
to overestimate polaron dissociation at the material interface. If, instead, model (C) is 
used, the obtained output current density characteristics is the dashed line in Fig. |9(a)] We 
observe a smoother trend than in previous cases for all applied voltages, and close to short 
circuit condition we note that the current density is further reduced since dissociation is 
assumed to occur only in the normal direction and on a significant portion of the interface 
E n is almost vanishing. In all the considered cases, the nonsmooth behavior at flat band 



condition (V apP i = 0.6 V) is to be ascribed to the discontinuity of d/3/dz at z = in ( 141. 

Fig. |9(b)| shows the results of the same analysis as above in the case of an exciton 
generation rate Q = 1.53 • 10 25 m _3 s _1 . The shape of the characteristics is very similar 
to those with low light up to a scaling factor of about 100, this suggesting a linearity 
between the output current density and the illumination intensity. Notice the absence 
of the bump for Vappi > 0.6 V in the case of model (B). This is a consequence of the 
increased magnitude of the charge carrier densities compared to the previously considered 
illumination that in turn determines stronger Coulomb attraction forces and hence more 
recombination phenomena. With reduced attractions, instead, charge carriers have more 
chances to escape from the interface following concentration gradients. 

In Fig. [10] we show the charge carrier densities in a device with geometrical data set to 
Lceii = 150 nm, L e z ec = 440 nm, Lr — 79 nm and Wr = 55 nm, at short circuit condition 
with exciton generation rate Q = 1.53- 10 25 m _3 s _1 . We first observe that computed 
charge carrier distributions in Fig. |10| 4eft) are in very good agreement with those of 
Fig. 3(i) in [43| and show the same peaks close to the vertical sides of the donor-acceptor 
interface. It is interesting to notice that the total number of holes in the donor material is 
higher than the number of electrons in the acceptor material because of the significantly 
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Parameter Symbol Numerical value 



Acceptor relative dielectric constant 




4 


Donor relative dielectric constant 


Sr,d 


4 


Temperature 


T 


298 K 


Poole-Frenkel mobility model parameters for electrons [4] 

Atn(O) 


3-10 m V s 
i.oo ■ iu v in 


Poole-Frenkel mobility model parameters for holes 


@ 

uJO) 

id 


1-10- 10 r^V^s" 1 
3- lCT 4 V- 1/2 m 1/2 


Exciton diffusion coefficient 


D e 


i • KnWs -1 


Exciton lifetime 


T e 


1 • 10 _9 s 


Exciton dissociation time 


Tdiss 


i ■ icr 12 s 


Polaron pair recombination rate constant 


krec 


1 ■ 10 8 s _1 


Singlet exciton recombination fraction 


V 


0.25 


Polaron pair zero-field dissociation rate constant 


kdiss (0) 


1 • 10 5 s _1 


Interface width 


2H 


2 • 10 _9 m 


Boundary condition parameters for electrons |17| 


a n 




1ms" 1 






P„ 


3.4995- 10 18 m- 2 s _1 


Boundary condition parameters for holes |17j 


k p 

Oip 




1ms" 1 






Pp 


3.4995- 10 18 m" 2 s" 1 



Table 2. Model parameter values used in the simulations of Sects. |5.2[ 

roiandroi 



different values of their respective mobilities. Negative charges can move through the 
device faster to be finally extracted at the cathode so that an overall positive charge 
builds-up in the device. The charge densities computed using models (B) and (C) exhibit 
a qualitatively similar profile with a gradual reduction of the magnitude compared to the 
result of model (A). This behavior is completely consistent with the previous analysis of 
the current-voltage characteristics predicted by the three models of kd ias . 

We conclude this preliminary validation of model (111 by illustrating in Fig. 11 the 
open circuit voltage V oc and short circuit current density J sc of a device with the same 
characteristics as in the previous set of simulations fo r value s of exciton generation rate in 
the range from 1.53 ■ 10 20 to 1.53 ■ 10 30 m -3 s - . Fig. 11(a) is in excellent agreement with 
Fig. 6(right) of [44], and indicates that models (A), (B) and (C) predict a linear behavior 
of V oc with respect to the logarithm of the exciton generation rate, as already pointed out 
in |43l 144] . Fig. 11(b) illustrates the current density J sc that can be extracted from the 
device at short circuit condition. The log-scale plot indicates that J 3C increases linearly 



in a wide range of illumination regimes until values of the order of 10 m s 



With 



more intense irradiation the increase becomes sublinear, suggesting that saturation of the 
device occurs due to more relevant excitonic and electron-hole recombination phenomena 
which in turn are a consequence of the increased densities. 
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Figure 10. Charge carrier densities [m 3 ] at short circuit condition 
with Q = 1.53- 10 25 m- 3 s _1 using models (A) (left), (B) (right) and 
(C) (bottom), respectively. 
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Figure 11. Open circuit voltage and short circuit current density as 
functions of the exciton generation rate. 
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5.3. The Role of Interface Morphology. In this section, we aim to investigate the 
role of interface configuration in affecting the OSC performance. Referring to Fig. [8] we 
set Lceii = L e i ec = 150 nm and Lr — 75 nm, and we analyze the importance of interfacial 
length by considering devices with an increasing density of interpenetrating structures, 
starting from a biplanar device and then taking decreasing values for the rod width Wr. 
Model parameters are the same as in the previous simulations and the exciton generation 
rate is Q = 1.53 • 10 25 m -3 s _1 . 




Interface Length [nm] 

Figure 12. Short circuit current density as a function of interface length. 

Fig. [12] illustrates the computed short circuit current density as a function of the inter- 
facial length for the various polaron pair dissociation rate models we previously considered 
in this section. In all cases, current saturation is predicted for high densities of nanostruc- 
tures due to the depletion of excitons in the interface area that in turn is a consequence of 
the abundance of dissociation sites. Computed saturation levels greatly differ among the 
three choices of the model for kdiss, in accordance with the analysis of Sect. |5.2| Fig. [12] 
also shows that when a biplanar device is considered, using model (C) a higher short 
circuit current density is obtained compared to the other approaches. An explanation 
of this result is that the electric field in this case is actually vertically directed and this 
fact, combined with the assumption that dissociation occurs only in the normal direction, 
brings to overestimate its rate (cf. the solid lines in Fig. [6p. Qualitatively similar results 
have been obtained in [8l l43l |44] , 

Also the orientation of the interface is expected to play a role in determining device 
operation and the following set of simulations aims to investigate this issue. This is a 
distinctive feature of our model that, to our knowledge, has not been treated in previous 
works. For a proper analysis, we allow the orientation of the donor-acceptor interface to 
change while its overall length remains almost constant, in order to single out the effect 
of the former and analyze it. 

The considered device geometry is shown in Fig. |13| where the number of rods is kept 
constant to four for each material but the incidence angle a is varied in a range from 90° to 
77° 11'. The geometric data are L ce u — L e i ec = 150 nm, Lr — 75 nm and Wr = 37.5 nm. 

Since the changes in the amplitude for a are small, the interface length does not 
vary significantly (as demonstrated by Fig. |14(a)[ | and we expect model (A) to be quite 
insensitive to such small modifications since E y mainly depends on the potential drop 
across the electrodes. Concerning with model (B), we again do not expect a relevant 
sensitivity to such variations in the interface morphology since the changes in E n and 
Et should balance in the overall contribution. We instead expect model (C) to be most 
sensitive since the normal field that is screened at the interface may experience significant 
variations as a function of the angle a. 
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Figure 13. Internal morphology with nanorods with a varying inci- 
dence angle a. 
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Figure 14. Interface length and short circuit current density as func- 
tions of a. 



Our expectations are confirmed by the results in Fig. 14(b) showing that the perfor- 
mance of the device in terms of computed short circuit current density does not vary with 
a when models (A) and (B) are considered, while if model (C) is used, an increase of the 
short circuit current density is observed as soon as the inclination of the nanorod structure 
is modified with respect to the initial configuration. This behavior can be explained as 



follows. The choice of model ( 15 \ predicts an increase of the dissociation for negative 
values of the normal electric field that is higher than the reduction for positive values of 
E n . Since at short circuit the electric field can be reasonably assumed to be directed along 
the y axis (i.e., from the cathode to the anode), the sides of each rod experience opposite 
normal fields. As a result, the overall effect is dominated by the contribution of the sides 
with negative fields and dissociation is enhanced. 

5.4. The Case of a Complex Interface Morphology. In this concluding section, we 
test the versatility of the model proposed in the present article in dealing with a very 
complex internal morphology as that shown in Fig. |15| In this regard, it is important 
to notice that the use of the microscale model |2|-(|5| would require an extremely fine 
grid resolution to accurately describe the volumetric terms in the active layer around 



the donor-acceptor interface, while the use of the macroscale model (111 has the twofold 
advantage of considerably simplifying the design of the computational mesh and reducing 
the size of the nonlinear algebraic system to be solved. 



Fig. 16(a) illustrates the computed charge carrier density at short circuit (V apP i — V) 



for the geometry of Fig. |15[ where the domain is a square 150 nm sided, with exciton 
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Figure 15. The computational mesh used to numerically solve the 
model of Sect. |3.4] in the case of a complex geometry. 
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Figure 16. Log-plot of charge carrier density [m ] and current-voltage 
characteristics for a device with very complex internal structure. 



generation rate Q = 1.53- 10 25 m~ 3 s _1 and using model (B) for kdiss- Notice in particular 
that the densities assume much higher magnitudes compared to those of Fig. [To] This is 
a consequence of the complexity of the geometry, where donor and acceptor form dead- 



end areas in which the charges are trapped and experience recombination. In Fig. 16(b) 
we show again a comparison of the current-voltage characteristics obtained using the 
three different polaron dissociation rate models. The differences among the obtained 
characteristic lines are reduced with respect to the previous simulated cases. In particular 
the computed short circuit current densities attain closer values with respect to more 
regular morphologies, such as that of Fig. [8j for comparable values of the interface length 
(approximately 900 nm), see Fig. 12 This is probably to be ascribed to the tortuosity of 



device internal morphology which makes interface recombination effects more significant 
than in the case of a more regular internal structure. 

6. Concluding Remarks and Future Perspectives 

The research activity object of the present article is a continuation of the mathematical 
study of organic photovoltaic devices started in [17] and is focused to: 

-: the accurate and computationally efficient modeling of photoconversion mecha- 
nisms occurring at the interface separating the acceptor and donor layers; 
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-: the investigation of the impact of the interface morphology and of polaron pair 
dissociation on device performance. 

With this aim, we propose a two-scale (micro- and macro-scale) multi-dimensional 
model for organic solar cell devices with arbitrary interface geometries. The microscale 
model is a system of incompletely parabolic nonlinear PDEs in drift-diffusion form set 
in a heterogeneous domain. The macroscale model is obtained via a micro-to-macro 
scale transition that consists of averaging the mass balance equations in the direction 
normal to the interface, giving rise to nonlinear transmission conditions parametrized by 
the interfacial width. This averaging procedure is similar to model-reduction techniques 
used in porous media with thin fractures [31], in reaction problems with moving reaction 
fronts 29 and in electrochemical transport across biological membranes [35] . The fact that 
in the macroscale model the interface is reduced to a zero width surface is further exploited 
to account for the local dependence of the polaron dissociation rate on the electric field 
orientation, which is the main advantage -together with the computational cost reduction- 
of our approach, as compared to previous multi-dimensional models [51 1441 l2l)] . 

Extensive numerical simulations of realistic device structures are carried out to study 
the performance of the proposed models and the impact of the lumping procedure. First, 
one-dimensional transient simulations under different working conditions are carried out 
to verify the accuracy of the macroscale model with respect to the microscale system. 
Results indicate that in the physically reasonable range of values for the parameter H the 
relative discrepancy between the micro and macroscale formulations is consistently below 
10%. Two-dimensional realistic device structures with various interface morphologies are 
then numerically investigated to assess the impact of our novel model for kdiss on the 
main device properties (short circuit current and open circuit voltage). Simulation results 
indicate that, if the electric field orientation relative to the interface is taken into due 
account, the device performance is determined not only by the total interface length but 
also by its shape. 



Research topics currently under scrutiny include: 



application of the proposed computational model to the study of more complex 
three-dimensional morphologies, as considered in [27] : 

investigation of more advanced models for carrier mobilities and polaron dissoci- 
ation rate, as well as the simulation of other material blends currently employed 
in the fabrication of up-to-date organic solar cells (see, e.g., [U0|9]); 
extension of the model to the general case where T + and T~ are free boundaries 
to be determined for each x G F and at each time level t > 0; 
a more thorough mathematical investigation of the proposed equation system (111 
in both stationary and time-dependent regimes to extend the analysis carried out 

in na. ' 
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Appendix A. Finite Element Discretization 

In this appendix, for sake of completeness, we provide more detail about the Finite Ele- 
ment (FE) discretization of the linear problem resulting from the application of time semi- 
discretization and linearization to the equation system as schematically described in 
Sectjl] Before we proceed we need to introduce some notation. 

The time semi-discretization consists of approximating the time derivative of the generic 
quantity U (representing any of the unknowns in ( |11[ )) as 

dU m=M 

(17) -7^- ~ WqUn + 2Z W m UN~m — W Un + d Nt M (U) , 



where ./V is the index of the current time step and M is the order of the adopted BDF 
formula. The notation djv,M (U) is used to group together the terms that depend only on 
results from past time steps and is therefore a known quantity at the JV-th time integration 
level. 

To treat the spatial discretization of the problem, we assume only for ease of pre- 



sentation that fi is a rectangular domain, as depicted in Fig. 3(b) but the approach 
remains completely valid also in the three-dimensional case, provided to replace "trian- 
gle" by "tethrahedron" and "edge" by "face". Let Th denote a conformal partition into 
open triangles K of the computational domain fi, h being the maximum diameter over 
all triangles, and let fi n ,ft and fi p ,h denote the finite element partitions of the subregions 
fi n and fi p , such that Fh — 9fi n ,h l~l 9fi p ,h is their separating interface, consisting of the 
union of a set of edges of Th- 

We introduce the following finite dimensional spaces of FE functions: 

(18a) V h = {v G C° (fi) , v\ K G ®x{K)VK G T h ) 

(18b) V h 9 = {v G V h , v = g at the nodes on T A U T c ,g 6 C° (T A U T c )} 
(18c) Vn,h = {v\ nnh ,veV h } 

(i8d) v p ,h = {v\ nph ,vev h } 

(18e) Vr,h = {v\ rh ,v€V h }. 

Let ipo G C° (r,4 UTc) be such that <pd\ Tc = and <pd\ Fa = V app i + Vu- Then, we 
denote by 

(19) y h = [e h , P h , n h , Ph ,ifh] T G (v£ x V r , h X V n , h X V p , h x V* D ) = Y h 

the vector of discrete unknown functions at a given quasi-Newton iteration of a given time 
step, and by 

(20) 5y h = [6e h ,5Ph,5n h ,5 Ph ,5ip h ] T G (Vh x V t ,h x V„, h x V p , h x T4°) = V h 

the corresponding increments to be computed in order to advance to the next iteration of 
the quasi-Newton method. 



The linear problem to be solved in order to compute the increments (201 reads: 
given y h G Y h , find Sy h G V h such that: 
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f D e VSe h -X7v+ f (—+wo]Se h v+ f \^— Se h — r/k r ec 5P h 
Jn Jn V r e / Jr h l T diss 



(21a) 



(21b) 

Li 



[ D e X7e h -X7v+ f f e (y h ) v + f g e (y h ) v 
Jn Jn Jr h 



2H 



Se h + (to + k di3S + k rec ) SP h - 2i?7 (p h Sn h + n h 8p h ) 



/ gpiyh) v, 

Jr h 



/ (D n , h V<5n h - n n ,h8n h V(ph) • Vv - / n n n h V8<p h • V« + / w Sn h v+ 
Jn n _ h Jn nh Jn n _ h 



/ \kdissSPh + 2H-/(p h Sn h + n h Sp h )\ v + / —8n h v 
Jr h L J Jr c K n 



(21c) 



- / (D„,/,Vn k - fJ,n,hn h V(ph) ■ V« + / fn(yh) V + / ffn(yh)u+ / fen(yfc) U 

J^n.h J^ n ,h Jr h Jr c 

/ (D p , h X78p h + fJ, p ,h8phVfh) -Vv + UpPhVSiph • V« + / too 5p fc u+ 

J^p.h J n p,h J n p,h 

/ [fc d i ss 5Ph + 2iJ7(p h< 5n h + n h< 5ph)l d + / — p h v = 
Jr h L J Jr A K v 

(21d) 

- / (Dp^Vrih + fi p ,hPh'V<Ph) ■ Vu + / fp{yh) v + / 5p(Yh)^+ / bp(yfe) w 



(21e) 



/ eVSifih -Vv + q I Sn h -q 5p h = - / eVip h • Vv + / yV(y/>) i> 
in J n nh J il p,h Lin in 



2(i 
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where v denotes in each of the equations (21 1 the test function in the appropriate FE 



space, and where 


we have made use of the followin 


i zza i 




= 1 — + wo Je h -Q + cLn,m (en ) 


(22b) 


9e (Yh) 


2H ^ p 


(22c) 


9p {yh) 


= e h + (WQ + k diss + k reC ) 


(22d) 


fn {yh) 


= w n h + d NiM (n h ) 


(22e) 




= kdissPh + 2Hjn h p h 


(22f) 


b n (yh) 


= — {a„n h — j3 n } 

fan 


(22g) 


f P (yh) 


= w p h + djv,A/ (ph) 


(22h) 


g P (yh) 


= kdissPh + 2Hjn h ph 


(22i) 


b P (yh) 


= — {otpPh - Pp} 

tip 


(22j) 


U (yh) 


= qn h - qp h . 



P h - 2H-fn h p h + dp 



Then, once system (211 is solved, the unknown vector is updated as 

yh ^yh + 8y h . 

A couple of further comments is in order about the numerically stable implementation 



of the FE linear system (211. First, note that Yi n ,h and fj,„,h in (21c I (respectively 
D p ,h and fi Pt h in ( 2 Id I ) are tensor diffusivities and mobilities chosen according to the 
Exponential Fitting stabilization technique as in [411 [3j 1181 1451 128] in order to avoid the 
onset of possible spurious oscillations in the discrete electron and hole densities due to 
drift terms. Second, all the integrals involving zeroth order terms are computed using 
the two-dimensional trapezoidal quadrature rule in order to end up with strictly positive 
diagonal (approximate) mass matrices 
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